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NONITERATIVE ESTIMATION OF A NONLINEAR PARAMETER* 

ARNE BERGSTROM 

ABSTRACT 

An algorithm is described which solves the parameters x = (aq, x 2 , ..., x m ) and 
p in an approximation problem Ax^y (p), where the parameter p occurs non- 
linearly in y. Instead of linearization methods, which require an approximate 
value of p to be supplied as a priori information, and which may lead to the find- 
ing of local minima, the proposed algorithm finds the global minimum by 
permitting the use of series expansions of arbitrary order, exploiting an a pri- 
ori knowledge that the addition of a particular function, corresponding to a 
new column in A, will not improve the goodness of the approximation. 


I. INTRODUCTION 

In the present paper the following approximation 
problem will be studied. Given 

(1) a function y(s, p) of an independent variable s (the 
extension to functions of several variables is trivial) and 
containing a parameter p which occurs nonlinear ly, and 

(2) a set of functions f^s), f 2 (s), .... f m (s) of s, 

determine a linear combination g(s, x) —x x f x {s) + x 2 f 2 {s) + 
... + x m f m (s) of the functions in (2), and a value of the 
parameter p in (1), such that the residual function 
e{s)=g(s, x)—y{s,p) over a certain interval in s is mini- 
mized in some sense. 

In discrete formalism, where the functions are given 
at discrete points s v s 2 , ..., s n , the approximation problem 
may for a large class of minimizing criteria be expressed 
as 

Ax=y(p)+e, (la) 

Be = 0, (lb) 


and B is a matrix which is derived from the criterion used 
for minimizing the residuals 



The most extensively used minimizing criterion is the 
least squares criterion, in which e T e is minimized. As 
is easily seen by differentiation, this criterion leads to 
a matrix B a equal to A T , the transposed coefficient 
matrix, with an additional row due to the nonlinear 
parameter and containing elements 




&y(8j, v ) 

dp 


Another example of a minimizing criterion is moment 
matching, where the first m + 1 moments of the right and 
left sides of Eqn. (la) are set equal. In the above formal- 
ism this is accomplished by using a matrix 


where A is a nxm matrix with elements a, 


x is the column vector 



B m = 


11 1 ... 1 
3 X 3 2 S 3 ... 3 n 

si si S 3 ... 4 
S? si si ... s\ 


\Si s 2 s 3 . 


y (p) is the column vector 


y(s 2 > P) 

\y{s n ,p)j 


The approximation problem may be illustrated by the 
following geometrical picture, see Fig. 1. We represent a 

* This work was in part performed at Southern Methodist Uni- 
versity, Dallas, Texas, on NASA Grant NsG 708. 
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Fig. 1. Geometrical representation of the nonlinear approxima- 
tion problem. 


function <p(s), given at the discrete points s lt s 2 , .... s n 

( <PM\ 

I as a point in a w-dimensional 

<p(s n )J 


space. Asp varies, the vector y (p) discussed above will in this 
representation describe a curve £(p); similarly, the left- 
hand side of Eqn. (la) will he represented by a 
m-dimensional hyperplane ~U A spanned by the columns 
of the matrix A, and the vector £ will connect a point 
on £(p) with a point on "U A . The criterion quantity 
which is used to measure the goodness of the approxi- 
mation, or in general how “close” two functions are to 
each other, defines in this picture a metric of the n- 
dimensional space. In the class of minimizing criteria 
expressed above as Eqn. (1) the metric is e r Me with 
A r M = B. This metric corresponds to the ordinary 
Euclidean metric e r € after an affine transformation 
characterized by a nxn transformation matrix N such 
that N r N =M, and is expressed in differential form by 
Eqn. (lb) with the implication that the shortest length 
of the residual vector € in the sense of the criterion used 
is when it is perpendicular to the hyperplane "U B which is 
spanned by the rows of B. 


II. PARAMETER SEARCH AND TAYLOR 
LINEARIZATION 

For the linear approximation problem of determining 
x for a fixed value of p from Eqn. (1), there exists an 
explicit and, if (BA) -1 is nonsingular, unique solution, 
x = (BA) -1 By. [To insure numerical stability, a somewhat 
different computational approach is, however, in general 
required (Householder, 1953; see also Golub, 1965).] 
The nonlinear part of the problem, i.e. the determination 
of the nonlinearly occurring parameter p, constitutes 
a somewhat more intricate problem both because there 
exists no such simple algorithm as in the linear case and 
because there may be several minima for the criterion 


quantity and ways must be found to ascertain that 
actually the lowest one is found. For this problem two 
classes of methods, employing parameter search or 
Taylor linearization, are usually used, either separately 
or in combination, and as a background for the following 
chapter a brief outline of these methods will here be 
given. 

In the parameter search, the linear approximation 
problem is solved for different fixed values of the non- 
linear parameter, and the corresponding measures of the 
goodness of the approximation are computed. In this 
way one aims at obtaining a coarse picture of the 
criterion quantity as a function of the nonlinear para- 
meter in order to discriminate between the region 
where the absolute minimum lies and regions where 
possible local minima may be found. In the former 
region a successively finer subdivision can then be made 
until the minimum is known within a required accuracy. 

In the linearization methods, a Taylor expansion of 
y(p) to the first order in p is performed around a value 
p 0 (which is assumed to be sufficiently close to the exact 
value), y(p) = y(p 0 ) + y'(p 0 ) A p + 0(Ap 2 ), where 


y'(Po) = 


dyfa. Po) 

dp 

dyfo, Po) 

dp 


dy{s n , Po) 

dp 


If we let A' denote a matrix which consists of the coeffi- 
cient matrix A as above with the additional column 


— y'(Po), a first approximation of the vector x' 



is obtained as x'= (BA') -1 By(p 0 ), where B = (A') r M. A 
new linearization and solution is then made around po 1 = 
p 0 +Ap (0> , and the procedure is repeated until A p lies 
within the required accuracy. 

Of the two classes of methods, the Taylor linearization 
is the fastest. A restriction is, however, that the initial 
guess must be sufficiently close to the exact value, 
otherwise there is always the risk that the method will 
find a local minimum. At the cost of a substantially more 
tedious computational procedure, the parameter search 
to a great extent excludes the possibility of finding local 
minima. However, if a too coarse scanning is made, there 
is also in this case a small risk that the absolute minimum 
may be lost. 

To illustrate the above difficulties we will study the 
following approximation problem. 


2 


FOA Reports, Vol. 7, No. 3, 1973 



Fig. 2. The function /(«) = 3 +a - exp (3s/2) + ®(0.9 sin (3«/2), 0.1) 
on the interval 0.6<«<1.0. Solid line represents least-squares 
fit. 


Example: Determine the parameters x lt x 2 and p in the 
expression 

g(s) =x 1 + x i s-e p, 

to give the best approximation in the least-squares sense 
to the points given in Fig. 2. (These points are computed 
as 

f(s) = 4 + 8 — e*' + <X»( — 0.9 sin 0.1), 

where 0(ra, a ) denotes a gaussian distribution with 
mean m and standard deviation a, which is included in 
order to simulate a disturbance with a systematic as well 
as a stochastic component.) 

Results: As is seen in Fig. 2, the function f(s) has a main 
trend which is convex from above. As long as p is 
sufficiently different from zero, this is also the case for 
the approximating function g{s), both if p is negative 
and if p is positive. When p becomes closer to zero, 
however, the approximating function degenerates to a 
linear function, and it is to be expected that the 
possibility to obtain a good fit to f(s) is less than with a 
convex function. This explains the appearance of a local 
minimum at p —2.5 in addition to the lowest minimum 
at p^\.5 in Fig. 3, where the criterion quantity Q 
(the sum of squares of the residuals) for the best choice of 
the corresponding linear parameter is plotted against the 
values of the nonlinear parameter p. This structure of the 
functional dependence of the criterion quantity may 
lead to difficulties when attempting to use the two algo- 
rithms discussed above to compute the nonlinear para- 
meter. If a negative value of p i3 used as starting value 
in the linearization method, the iterative procedure will 
in general find the local minimum instead. This turns 
out to be true also for more elaborate linearization 
methods such as those proposed by Davidon (1959) and 
Powell (1965). Due to the sharpness of the lowest mini- 
mum there may be some danger that also the parameter 



Fig. 3. Sum of squares of residuals Q versus the value of the non- 
linear parameter p. (Note the suppressed zero.) 

search misses this minimum if a too coarse initial scanning 
is performed. 

In the new, noniterative algorithm which will be 
discussed in the following, these difficulties are removed. 

III. THE NEW ALGORITHM 

In contrast to the Taylor linearization method, which 
requires an approximate value of the nonlinear para- 
meter as a priori information, the new algorithm which 
we will develop here exploits a possible a priori knowledge 
that the addition of a particular function to the 
left-hand side of Eqn. (la) will not improve the goodness 
of the approximation, i.e. that the vector in the above 
geometrical model which corresponds to the additional 
function is orthogonal, in the sense of the criterion, to 
the residual vector e. Especially if the dimensionality 
of the space is high, as is desired for reasons of numerical 
accuracy, there is a multitude of functions obeying this 
orthogonality requirement and, as we shall see in 
Chapter IV, experience also shows that there is as a 
rule no difficulty in finding such functions; suitable 
functions are even in many cases suggested by the very 
assumptions made in formulating the problem, e.g. that 
certain terms can be neglected in the mathematical 
description of the actual problem. 

Using the picture of Chapter I, the original approxima- 
tion problem corresponds to the geometrical problem of 
finding a point on the curve C(p) and a point on the 
hyperplane "U A such that the residual vector e, orthogonal 
to 7t B , is minimized. By including a new function (as 
column ra+ 1) in A as discussed above, a new hyperplane 
"U'a’ now °f dimension m + 1, is formed (Fig. 4). Due to the 
assumptions inherent in the choice of the additional 
function, this hyperplane has the property of being 
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Fig. 4. The nonlinear parameter is solved from the best intersec- 
tion between the hyperplane (in the figure one -dimensional) 
and the projection of c(p) on the hyperplane 


orthogonal in the sense of the metric to the residual 
vector c above. 

The solution to Eqn. (1), x = (BA) -1 By(p), for brevity 
called c (p) in the following, gives x as a function of p 
corresponding to a curve on the hyperplane fl' A which is 
the projection of C(p). (We here assume that the vector 
corresponding to the additional function is not ortho- 
gonal to C(p)-) As follows immediately from the above 
constructions, one of the intersections between the 
projected curve and the original hyperplane "U A is 
identical with the endpoint of the residual vector e, 
which we wanted to determine. This implies that the 
value of the nonlinear parameter can be solved by 
setting the coefficient x m+l for the additional function 
equal to zero, i.e. the nonlinear parameter is one of the 
roots to the equation c m+1 (p) =0. This expression is 
greatly simplified if y (p) is given in the form of a series 
expansion. 

If the functions y (p) are given as power series expan- 
sions, assuming the power series expansion of c m+1 (p) to be 
convergent, the new algorithm can be expressed in 
algorithms for the solution of linear approximation 
problems and polynomial roots by using the following 
lemma. 

Lemma. Let a linear approximation problem with n 
observations and m unknown parameters be given as 

Ax = y+£, 

Be = 0, 

where (BA) -1 is nonsingular, so that there exists a unique 
solution. If we perform an arbitrary decomposition of 
y into y=yi+y a + ...+y„ the solution x to the original 
problem may be written as x=x 1 +x 2 + ...+x,, where 

4 



Fig. 5. Geometrical interpretation of the decomposition lemma. 


x„ t = l, 2, .... I, are the solutions x, = (BA) -1 By, to 
the decomposed approximation problems 

Ax, = y,+e„ 

Be, = 0. 

Proof. In each of the decomposed approximation 
problems, Ax,=y, +e„ Be,=0, we have n+rn equations 
which, since (BA) -1 is assumed to be nonsingular, com- 
pletely determine the n + m unknowns x, and e,. By 
summing over i, i — 1, 2, ... I, we obtain 

A(x 1 +x 2 +...+x,)=y 1 +y 2 + ...+y,+€ 1 +e 2 +...+e„ 
B(e 1 +e 2 + ... +e,) = 0, 

and, since yi+y 2 + ...+yi=y, we obtain 
x = x 1 +x 2 + ...+x,. Q.E.D. 

In the geometrical picture described in Chapter I this 
lemma has a simple interpretation, see Fig. 5. The 
vector Ax is the projection of the vector y on the 
hyperplane "U A . The lemma merely states that when 
y is expressed as a sum of vectors y x , y 2 , .... y ( , the vector 
Ax may be expressed as the sum of the projections of 
yi, y 2) — > Ji on -Ha, i-e. the vectors Ax,, Ax 2 , ..., Ax,. 

Returning now to the nonlinear approximation prob- 
lem discussed earlier, 

Ax=y (p)+e, 

Be =0, 

where the coefficient for the artificial dimension in- 
troduced is x m+1 , we write the first system of equations 
in explicit form 

a ii x i 'h%2 a '2“ r = yio^ViiP Jr yi2'P i + ••• +£i 

a 21 X 1 + a 22 X 2 + . . . + «2 . CT+1 *m+ 1 = Z /20 "ky22jP 2 + ~i" £ 2 

a nl x l ~l~ a n2 x 2~t~ ••• + a n,m+l x m+l ~ V V n\V ^ V n^V 2, ^ ■ 

+ £„. 

In correspondence with the power series expansion in the 
FOA Reports, Vol. 7, No. 3, 1973 



right-hand side, we now decompose x v x 2 , x 3 , ..., x m+l 
as x 1 =x 10 + x n + x 12 + x 2 =x 20 +x 21 +x 22 + ■■■, etc. After 
introduction of the new unknowns x[o = x 10 , x'n —x^p, 
*12 x^ 2 lp 2 t •••» *20 — x 2 q , x 2 i~x 21 Ip, x 2 2~x 22 jp 2 i etc., 
we may from the overdetermined equation system above 
form a new set of overdetermined equation systems: 

^11 *10 +" ®12 *20 + . . . + oq.m+i * m+1,0 = 2/l0 P" e 10 
a 21 *10 + a 22 *20 + . . . + a 2 .m+l * m+1,0 = 2/20 d" £20 

a nl *10 + a n2 *20 + • • • + ® n ,m+l *m + l,0 = VnO d" £ n 0> 

from which the unknowns xj = (*io,*2o> ■•••*m+i.o) can be 
determined; 


An alternative to the above geometrical approach of 
describing the algorithm is as follows. 

The a priori assumption is equivalent to saying that 
the approximation problem 

(A a m+ j) r )=y(p)d-E, 

\*m+l/ 



has the solution x m+1 = 0. Eliminating e we obtain 



Ba, 

b r 


'm+1 

m+l a m+l. 




<hi *u p d- <h 2 *21 p +... + a 1 _ m+1 x' m+lil p = y n p + £ U 
«21 *11 P + «22 *21 P + ■ ■ ■ + «2.m+l *m+l,l P = VilP+ £21 

«7il*nP + ®,.2*2lP + . . . + a n .m+l *m+l,l p — y nl p+ e nl> 

from which, after division by p, the unknowns xj = [x n , 
*2i> •••> *m+i , 1) can be determined (since all residuals 
£ u> e 2i> •••> £ ni are affected equally by the divison with p, 
the division does not change the problem); 

a n x ' 12 p 2 + a l2 x 22 p 2 + ... + a lmn x ' m+li2 p 2 = y 12 p 2 + e 12 

«2i z'12 P a d- a 22 x^ p 2 + . . . + « 2 , m+1 *^ +12 p 2 = y 22 p 2 + e 22 

a„i x'i 2 p 2 + a n 2 x 22 p 2 + ...+ a nm+1 x ' m+12 p 2 = y n2 p 2 + e n2 , 

from which after division by p 2 the unknowns x 2 = 
(*12, *22) -i *m+i.2) can be determined; 

etc. for the other powers of p. 

Since only the right-hand side differs between the 
cases, the solution of these equation systems is a very 
fast procedure. 

According to the decomposition lemma the solution 
vector x to the original problem is given as 

x^Xo + Xi + Xa-t- ... =xi + xjp+x 2 p a + — 

Following the discussion in the beginning of this chapter, 
the value of the nonlinear parameter can then be solved 
from the equation a; m+1 =0, selecting the best root in 
the sense of the criterion used. Once the value of p is 
known, x 1 , x 2 , ..., x m can also be computed as x 1 = 
Xi 0 +x'xiP+x'i2P 2 + * 2 = *20+*2lP+*22P 2 + ---i etc. 

To illustrate the conciseness with which the new 
algorithm can be programmed using existing algorithms 
for linear least-squares problems and polynomial roots, 
and to provide a model for its implementation in other 
programming languages, a FORTRAN program of a least- 
squares version expressed in subroutines from the IBM 
Scientific Subroutine Package (IBM, 1970) is given in 
Fig. 6. 


Assuming the matrix (BA) to be nonsingular, x m+1 can 
be solved from this equation, e.g. by the method of 
bordering (Faddeeva, 1959), i.e. 

*m+l b m +i(I P a / b ) y(P)t *771 b m + j(I - P^/s) a m+l» 

a m 

where 

P 4/a = A(BA)- 1 B 

is a projection operator, and we have assumed that 
a m + 0. If we define r m+1 by 

r m+i = (I — h m+1 

(note that r m+1 is easily computed as the residual vector 
corresponding to the approximation problem with right- 
hand side b m+1 ), we get the simple formula 

*771+1 “ r m+iy(p)l *771 = r 771+l a 771+l- 

*771 

The condition * m+1 =0 now gives a nonlinear equation 
from which p can be solved. In particular, if y(p) is given 
by a power series expansion y(p) = y 0 + yiP +y2P 2 + ••■» 
we get a 0 +a 1 p+a 2 p 2 + ...= 0 , where a k =vl + 1 y k . 

In the case of the least-squares criterion, P^ /B is sym- 
metric, and the formula for r ra+1 reduces to 

r i7i+i = P — A(A r A) -1 A r ] b m+1 . 

As above, r OT+1 is here easily computed (as the residual 
vector corresponding to the least-squares problem with 
right-hand side b m+1 ). 

The method can be characterized by saying that an a 
priori vector a m+1 is used instead of — dy/dp in the (m+ 1)- 
th column in A. To guard against an unfortunate choice of 
the artificial dimension, one can after the first value p =p* 
has been obtained choose a m+1 = — {dyldp) v ~ D * as new vec- 
tor and recompute p. This method can of course easily 
be iterated, and since only the right-hand side changes, 
it is a very fast procedure. 

We will here also give the expressions for the statistical 
uncertainties in the estimates of the parameters under 
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SUBROUTINE TGTHPL 
PURPOSE 

LEAST SQUARES SOLUTION OF THE VECTOR X AND THE PARAMETER P 
IN AN OVERDETERMINED EQUATION SYSTEM WRITTEN IN MATRIX FORM 
AS A*X=Y(P), WHERE P OCCURS NON-L INEARLY IN THE VECTOR YIP) 

USAGE 

CALL TGTHPL! A , V,S ,X , R , AX, Y ,Q ,E, W1 , W2 , N ,M, L , I ER1 , 1 ER2, I ER3 ) 
DESCRIPTION OF PARAMETERS 

A - THE N BY M INPUT MATRIX OF COEFFICIENTS FOR THE 

LINEAR UNKNOWNS,. THE M-TH COLUMN REPRESENTS THE 
ARTIFICIAL DIMENSION USED TO FORM THE TANGENT 
HYPERPLANE, SEE UNDER METHOD BELOW 

V - THE N BY L INPUT MATRIX OF COEFFICIENTS IN THE 

•TAYLOR EXPANSION OF THE RIGHT HAND SIDE IN POWERS 
OF THE .NON-LINEAR UNKNOWN, ORDERED FROM LCW TO 
HIGH ORDER 

S - OUTPUT VARIABLE CONTAINING THE SUM OF SQUARES OF 
THE RESIDUALS 

X - OUTPUT VECTOR OF LENGTH M CONTAINING THE SOLUTION. 

LOCATION M CONTAINS THE NON-LINEAR UNKNOWN 
R - CONTAINS ON RETURN TEE M BY L MATRIX OF COEFFI- 

CIENTS IN THE EXPANSION OF THE LINEAR UNKNOWNS 
AX - WORKING STORAGE OF LENGTH MAXI2*N,U. THE FIRST 

N LOCATIONS CONTAINS ON RETURN THE LEFT HAND SIDE 
VECTOR A*X 

Y - CONTAINS ON RETURN TEE RIGHT HAND SIDE VECTOR Y 

OF LENGTH N 

Q - VECTOR OF LENGTH L-l CONTAINING ON RETURN THE REAL 

PARTS OF THE ROOTS TO THE POLYNOMIAL EQUATION FOR 
THE NON-LINEAR UNKNOWN 

E - VECTOR OF LENGTH L-l CONTAINING ON RETURN THE 

IMAGINARY PARTS OF THE ROOTS TO THE POLYNOMIAL 
EQUATION FOR THE NON-LINEAR UNKNOWN 
W1 - WORKING STORAGE OF LENGTH KAX<N*M,L> 

W2 - WORKING STORAGE OF LENGTH N*L 

N - NUMBER OF ROWS IN EQUATION SYSTEM 

M - NUMBER OF COLUMNS IN EQUATION SYSTEM I INCLUDING 

ARTIFICIAL DIMENSION) 

L - NUMBER OF TERMS IN TAYLOR EXPANSION OF YIP) 

I ERl - ERROR MESSAGE FROM LINEAR LEAST SQUARES SUBROUTINE 
I ER2 - ERROR MESSAGE FROM POLYNOMIAL ROOTS SUBROUTINE 
I ER3 - ERROR MESSAGE FROM TGTHPL 
I ER3=0 - NO ERROR 

I ER3= 1 - NO REAL ROOTS. THE RESULTS GIVEN IN 

THIS CASE ARE COMPUTED USING THE BEST 
REAL PART OF THE ROOTS 

REMARKS 

MATRICES A , V, R ARE GENERAL MATRICES 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 

IBM SYSTEM/360 SCIENTIFIC SUBROUTINE PACKAGE ROUTINES: 
MCPY,RCPY,GMSUB,GMPRD,MATA,LLSQ,PVAL,POLRT 
(REFERENCE: IBM PUBLICATION C-H20-0205-4 ) 

METHOD 

THE ROUTINE USES THE NON-ITERATIVE ALGORITHM DESCRIBED BY 
A. BERGSTROM, FOA REPORT B4059-M4 (1973). 

A VECTOR CORRESPONDING TO A FUNCTION WHICH WILL NOT IMPROVE 
THE GOODNESS OF THE APPROXIMATION HAS TO BE SUPPLIED BY 
THE USER AS COLUMN M IN MATRIX A. THIS VECTOR IS USED 
TO FORM A TANGENT HYPERPLANE UPON WHICH THE CURVE Y(P) IS 
PROJECTED. THE SOLUTION TO THE ORIGINAL PROBLEM IS THEN 
OBTAINED BY SETTING THE COEFFICIENT FOR THE ADDITIONAL 
VECTOR EQUAL TO ZERO 
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Fig. 6. FORTRAN IV program print-out of least-squares version of the algorithm. 
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SUBROUTINE TGTHPL < A , V ,S , X, R , AX, Y ,Q ,E , WI , W2 , N, M, L, I ERI , I ER2, I ER3 > 
DIMENSION AI 1I,V(1 J ,X(I), R(1),AX( 1 1,Yll) ,QC 1) ,E(1J tWlUI, W2U» 
CALL MCPY(A,W1,N,M,0) 

CALL MCPY(V,W2,N,L,0) 

CALL LLSQIW1,W2,N,M,L,R,Y,0. , IER1, AX > 

CALL RCPY(R,M,W1,M»L»0) 

IR=L-1 

CALL POLRT(Wl,W2, IR,Q,E, IER2) 

I ER3= 1 

DO 10 1=1, IR 

CALL EVALCI , Q,R, X, V ,Y,A , AX ,S ,N, M,L ,W1,W2J 
IFI I.NE. l.AND.t IER3«NE. loANb.tEI I) 0 NE 0 O 0 0 OR 0 S 0 GE 0 SO )«0Ro 
& E ( I ) .NE.O. . AND.S.GE, SO) ) GOTO 10 
IFIE(I)* EQ. 0. ) I ER3=0 
10=1 
SO=S 

10 CONTINUE 

CALL EVAL(IO,Q,R,X,V,Y, A, AX , S, N , M, L, Wl, W2J 

X ( M ) = Q1 1 0 J 

RETURN 

END 

SUBROUTINE EVAL ( I , Q , R ,X , V , Y, A, AX , S ,N , M, L,W1 , W2 » 

DIMENSION Q(l),R(l),Xm,vm,Ym,A(l),AX(l),Wl(l),W2(l> 

DO 10 J=1,M 

CALL RCPYIR, J,W1,M,L,0> 

10 CALL PVAL(X(J),Q<I),W1,L) 

DO 20 J=1»N 

CALL RCPYIV, J,W2,N,L,0J 
20 CALL PVALIYI J>,Q1I) ,W2,L) 

CALL GMPRDI A,X, AX,N,M,11 
CALL GMSUBI Y,AX,W1,N»1» • 

CALL MATA(W1,S,N,1,0) 

RETURN 

END 


TGTH 750 
TGTH 760 
TGTH 770 
TGTH 780 
TGTH 790 
TGTH 800 
TGTH 810 
TGTH 811 
TGTH 820 
TGTH 830 
TGTH 840 
TGTH 850 
TGTH 660 
TGTH 870 
TGTH 880 
TGTH 890 
TGTH 900 
TGTH 910 
TGTH 920 
TGTH 930 
TGTH 940 
TGTH 950 
TGTH 960 
TGTH 970 
TGTH 980 
TGTH 990 
TGTH1000 
T GTH1010 
TGTH1020 
TGTH1030 
TGTH1040 
TGTH1050 
T GTH1060 
TGTH1070 


Fig. 6 (continued). 


the assumption that the e, are independent, have zero 
mean, and standard deviation a. As discussed above, the 
solution to the approximation problem is given as 
x = (BA) _1 By(p), with * m+1 = 0. Denoting a particular 
estimate by x* and the matrix (BA) -*8 by D, we im- 
mediately obtain x* = D(Ax + e)=x + De. By definition, 
the variance-covariance matrix of x* is 

V(x*) = E[(x* -x) (x* — x) T ], 
i.e. 

V(x*) = <r 2 DD r . 

In the least-squares case B=A r , and V(x*) reduces to 
the familiar expression V(x*) =a 2 (A r A) _1 . Since x m+1 is 
forcibly set to zero, an uncertainty Ax* m+1 in x* m+l 
corresponds directly to a change — Ax* m+1 in the 
constant term in the Taylor expansion of x m+1 (p). Such 
a change in the constant term is related to a change in 
p by an amount A p from the computed value p 0 , given as 

A p~ ^* +1 . 

x m+i(Po) 

IV. APPLICATIONS 

In this chapter we will discuss some applications and 
limitations of the proposed algorithm. 


(1) When using the program given in Fig. 6 on the 
example in Chapter II, which as we saw proved to be 
a rather difficult task for the customary methods, no 
difficulties were encountered at all, and the program 
gave a very fast and accurate determination of the 
lowest minimum. If, for instance, a ninth order Taylor 
expansion around p= 0 was performed and the function 
s 2 was used to form the artificial dimension, the value 
obtained for the nonlinear parameter differed only by 
about 0.1 % from the exact value of the lowest minimum, 
a difference which of course is far within the uncertainty 
of the order of 5 % in p inherent in the problem due to 
the statistical scatter in the input data. To illustrate the 
statement in Chapter III that there is a multitude of 
possible functions available for the construction of the 
artificial dimension, parameter estimates were also made 

i 

using the functions s', 1/s 1 , j/s^ cos's, In' s, with 1 = 3, 
4, 5, 6, 7 to form the artificial dimension; in all 
these cases the value obtained for the nonlinear para- 
meter was far within the statistical fluctuations in p 
just as in the case of s 2 above. 

(2) The algorithm assumes, in the formulation of Chapter 
III, the functions y (p) to be given as Taylor expansions 
of arbitrary order in p, i.e. y{p) = y 0 + yiP + y 2 P 8 + — - 
In some cases an expansion of the type 
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y(j , )“p+^Fx+ •••+yj+yi+i?’+y!+22 ,a + ••• 

may be more suitable to approximate the functional be- 
havior of y (p). This case is easily reduced to the one 
previously discussed by multiplying both sides in the 
approximation problem, Eqn. (1), by the unknown 
constant p l (since p l is constant, this merely corresponds 
to a change of all weights by the same amount), after 
which the algorithm can be used to calculate p l \ and p. 

(3) An important limitation to the applicability of the 
algorithm in the form presented here is the requirement 
that the problem can be brought on the form Ax = 
y(p)+c, i.e. to separate the linear and nonlinear para- 
meters. This means that an approximation problem of the 
type, say, f(s)^a + bs+ce xs +de xs \ where we want to 
determine a, b, c, d, and a, cannot be solved using the 
algorithm described here. (A problem of the type g(s) 

a + bs + ce xs is, however, after division by c easily 
reduced to a form g{s)fc — a/c — bsjc^e 13 , for which the 
algorithm is applicable.) 

(4) Potential limitations to the applicability of the algo- 
rithm are also inherent in the construction of the artificial 
dimension as discussed in Chapter III and the con- 
vergence properties of the power series expansion of the 
right-hand side. For reasons of numerical accuracy 
alone, one aims in general to have a number of observa- 
tions which is large compared to the number of para- 
meters to be determined; to find a function which does 
not improve the approximation will then as a rule, and as 
we also saw in paragraph (1) above, offer no problem. 
Also the convergence properties of the series expansion 
is in practice seldom any serious limitation, since the 
algorithm uses the projection c(p) of y (p) on “U' A rather 
than y(p) itself, and the projection has in general a more 
tranquil behavior than the original function. 


V. CONCLUDING REMARKS 

The noniterative algorithm presented in this paper 
might, if used within the limitations discussed in 
Chapter IV, be a useful alternative to the customary 
methods for solving approximation problems containing 
one nonlinear parameter. 

In a general problem with several linear and nonlinear 
parameters 


A(p)x = y(p)+e, 

B(p)e = 0, 

the solution of the linear parameters x should always in 
virtue of their uniqueness be separated from the solution 
of the nonlinear parameters p. Especially simple is this in 
the case of gradient-free algorithms such as the one 
given by Powell (loc. cit.). In this case the criterion 
quantity 

Q(x, p) = [A(p)x-y(p)fM[A(p)x-y(p)] 

can be reformulated as 

<2(P) = [A(p)x*(p) — y(p)] r M[A(p)x*(p) — y(p)], 

where the solution of the linear parameters, 

x*(p) = [B(p)A(p)]- 1 B(p)y(p), 

is calculated in the function-evaluating routine in the 
program. Here the present algorithm might be introduced 
instead to calculate both x* and one of the nonlinear 
parameters, thus reducing the dimensionality of the 
nonlinearity, which may lead to a substantial increase in 
safety as well as computational economy. 
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